中子飞行时间谱测量原理¶
探测器参数¶
探测器的几何设置见下图：
- 探测器长度：2L, 忽略探测器宽度。
- 探测器厚度：$\Delta D$
- 光在闪烁体的传播速度: $v_{sc}$
- 光在闪烁体的衰减常数：$\lambda$
- 粒子从起始位置到探测器中心的最小距离：D
- 粒子入射位置：x
- 飞行距离：$d=\sqrt{D^2+x^2}$
飞行时间：以靶为起点
中子飞行时间：$TOF_n(ns)=\frac{72 \cdot d(m)}{\sqrt{E_n(MeV)}}$
$\gamma$ 飞行时间：$TOF_{\gamma}(ns)=3.333\cdot d(m)$
探测器两端时间信号(不考虑其他offset)：
- $t_u=TOF+(L-x)/v_{sc}$
- $t_d=TOF+(L+x)/v_{sc}$
- $x=(t_d-t_u)*v_{sc}/2$
- $TOF=(t_u+t_d)/2 -L/v_{sc}$
探测器两端能量(x 处能量沉积$Q_0$):
- $Q_u=Q_0 exp[-(L-x)/\lambda]$
- $Q_d=Q_0 exp[-(L+x)/\lambda]$
- $x=\frac{\lambda}{2}ln\frac{q_u}{q_d}$
- $Q_0=e^{2L/\lambda}*Q_u*Q_d$
代码¶
// 运行方法：
// 将下列代码保存到tree.cc
// 在ROOT环境中运行：
// .x tree.cc
// 生成 tree.root 文件
void tree(){
//常量声明
const Double_t D=500.;//cm, distance between target and the scin.(Center)
const Double_t L=100.;//cm, half length of the scin.
const Double_t dD=5.;//cm, thickness of the scin.
const Double_t TRes=1.;//ns, time resolution(FWHM) of the scintillator.
const Double_t Lambda=380.;//cm, attenuation lenght of the scin.
const Double_t QRes=0.1;//relative energy resolution(FWHM) of the scin.
const Double_t Vsc=7.5;//ns/cm, speed of light in the scin.
const Double_t En0=100;//MeV, average neutron energy
const Double_t EnRes=50.;//MeV, energy spread of neutron(FWHM)
const Double_t Eg0=1;//MeV, gamma energy
const Double_t Rg=0.3;//ratio of gamma,ratio of neutron 1-Rg
//1. 定义新ROOT文件，声明新的Tree
TFile *opf=new TFile("tree.root","recreate");//新文件tree.root，指针 *opf
TTree *opt=new TTree("tree","tree structure");//新tree，指针 *opt
//2. 声明在tree结构中定义需要的变量分支
Double_t x;//入射位置
Double_t e;//能量
int pid; //粒子种类，n:pid=1,g:pid=0
Double_t tof, ctof;//TOF:粒子实际飞行时间，cTOF：计算得到的TOF
Double_t tu, td;
Double_t qu, qd;
Double_t tu_off=5.5;//time offset，//PMT的度越时间+电缆上的传输时间
Double_t td_off=20.4;//time offset
//3. 将变量地址添加到tree结构中
//第一个参数为变量名称，第二个为上面定义的变量地址，第三个为变量的类型说明，D表示Double_t。
opt->Branch("x", &x, "x/D");
opt->Branch("e", &e, "e/D");
opt->Branch("tof", &tof, "tof/D");
opt->Branch("ctof",&ctof,"ctof/D");
opt->Branch("pid", &pid, "pid/I");
opt->Branch("tu", &tu, "tu/D");
opt->Branch("td", &td, "td/D");
opt->Branch("qu", &qu, "qu/D");
opt->Branch("qd", &qd, "qd/D");
// histogram，ROOT文件中除了TTree结构外，还可存储histogram，graph等
TH1D *hctof=new TH1D("hctof","neutron time of flight",1000,0,100);
TRandom3 *gr=new TRandom3(0);//声明随机数
//4. 循环，计算变量的值，逐事件往tree结构添加变量值。
for(int i=0;i<100000;i++){
x=gr->Uniform(-L, L);//均匀入射在中子探测器表面.
Double_t Dr=D+gr->Uniform(-0.5,0.5)*dD;//粒子在探测器厚度范围内均匀产生光信号
Double_t d=TMath::Sqrt(Dr*Dr+x*x);//m, flight path
if(gr->Uniform() < Rg) { //判断为gamma入射
pid=0;
e=Eg0;
tof=3.333*(d*0.01);
}
else { //neutron
pid=1;
e=gr->Gaus(En0, EnRes/2.35); // neutron
tof=72./TMath::Sqrt(e)*(d*0.01);//ns
}
tu=tof+(L-x)/Vsc+gr->Gaus(0,TRes/2.35)+tu_off;
td=tof+(L+x)/Vsc+gr->Gaus(0,TRes/2.35)+td_off;
ctof=(tu+td)/2.;//simplified calculation.
hctof->Fill(ctof);
Double_t q0=e*gr->Uniform();//energy of recoil proton in plas. 0-En
qu=q0*TMath::Exp(-(L-x)/Lambda);
qu=gr->Gaus(qu,qu*QRes/2.35);
qd=q0*TMath::Exp(-(L+x)/Lambda);
qd=gr->Gaus(qd,qd*QRes/2.35);
opt->Fill();//5.将计算好的变量值填到Tree中
}
// 6.将数据写入root文件中
hctof->Write();
opt->Write();
opf->Close();
}
练习：¶
- 理解上述1.-6.步骤的逻辑关系。自行练习如何将数据已TTree格式存到ROOT文件中。
- 理解如何在模拟中加入探测器分辨率信息。
查看ROOT文件¶
- juypter运行tree->Draw()时，需要有网络链接才能正常显示。
In [1]:
%jsroot on
In [2]:
TFile *ipf=new TFile("tree.root");//root -l tree.root
ipf->ls()//ROOT 环境下 > .ls
In [3]:
//观察tree的结构
tree->Print()
观察指定事件数据¶
In [15]:
tree->Show(0);//显示指定的事件数据 event=1
In [5]:
tree->Show(100);
In [6]:
TCanvas *c1=new TCanvas();//* 在ROOT环境下可省略
c1->Clear();//* 在ROOT环境下可省略
tree->Draw("td-tu>>htx(500,-20,50)");//位置一维分布
c1->Draw();//* 在ROOT环境下可省略
In [7]:
tree->Draw("tof>>hh(1000,0,100)");//实际飞行时间
c1->SetLogy();
c1->Draw();
In [8]:
TH1D *hh=(TH1D*) ipf->Get("hctof"); //在代码中得到文件内数据指针的方法，在ROOT环境下课直接用 hctof->Draw()
hh->Draw();//显示文件内histogram
c1->Draw();
思考题1：观察上两个图有什么不同，分析其原因¶
In [13]:
c1->SetLogy(0);
gStyle->SetPalette(1);
tree->Draw("ctof:td-tu>>(2000,-20,50,1000,0,100)","","colz");//二维图显示
c1->Draw();
检查参数之间的关系是否与物理条件一致¶
- gamma的到达探测器的时间(飞行时间)与探测器的x位置无关？
In [14]:
tree->Draw("ctof:td-tu>>(2000,-20,50,50,40,45)","pid==0","colz");//加入条件,gamma
c1->Draw();//为什么在两侧上弯？
In [11]:
tree->Draw("td-tu:x","","colz");//观察两个参量的关联
c1->Draw();
In [12]:
tree->Draw("log(qu/qd):x","","colz");
c1->Draw();
In [3]:
!jupyter nbconvert 1.1_create_tree.ipynb --to html
In [ ]: